;simply traces 4 axially alligned N pixel length rays and estimates crown radius based on measurements
;Roberts et al (2005)
;if crown height on a ray goes below seed height - hthresh the ray ends
;return is a list of radii associated with input seed points
;seeds is 2 by M array
function traceraycircle4, chm, n, seeds, hthresh, rays=rays
	numseeds = n_elements(seeds[0,*])
	radii = ulonarr(numseeds)
	length = intarr(numseeds) ;will be averaged
	if (arg_present(rays)) then r = dblarr(n,numseeds, 4)

	;up
	arrin = ulonarr(2,n,n_elements(seeds[0,*]))
	arrin[0,*,*] = rebin(seeds[0,*],n,numseeds)
	arrin[1,*,*] = rebin(seeds[1,*],n,numseeds) + rebin(lindgen(n),n,n_elements(seeds[0,*]))

	for i=0L, numseeds - 1 do begin
		chmA = chm
		chmA[where(chm lt chm[seeds[0,*], seeds[1,*]]) - hthresh]=max(chm)
		val = min(chmA[arrin[0,*,i],arrin[1,*,i]], minsc)
		length[i] += minsc
		if (arg_present(rays)) then r[*,i,0] = chm[arrin[0,*,i], arrin[1,*,i]]
	endfor

	;down
	arrin = ulonarr(2,n,n_elements(seeds[0,*]))
	arrin[0,*,*] = rebin(seeds[0,*],n,numseeds)
	arrin[1,*,*] = rebin(seeds[1,*],n,numseeds) - rebin(lindgen(n),n,n_elements(seeds[0,*]))

	for i=0L, numseeds - 1 do begin
		chmA = chm
		chmA[where(chm lt chm[seeds[0,*], seeds[1,*]]) - hthresh]=max(chm)
		val = min(chm[arrin[0,*,i],arrin[1,*,i]], minsc)
		length[i] += minsc
		if (arg_present(rays)) then r[*,i,1] = chm[arrin[0,*,i], arrin[1,*,i]]
	endfor

	;left
	arrin = ulonarr(2,n,n_elements(seeds[0,*]))
	arrin[1,*,*] = rebin(seeds[1,*],n,numseeds)
	arrin[0,*,*] = rebin(seeds[0,*],n,numseeds) - rebin(lindgen(n),n,n_elements(seeds[0,*]))

	for i=0L, numseeds - 1 do begin
		chmA = chm
		chmA[where(chm lt chm[seeds[0,*], seeds[1,*]]) - hthresh]=max(chm)
		val = min(chm[arrin[0,*,i],arrin[1,*,i]], minsc)
		length[i] += minsc
		if (arg_present(rays)) then r[*,i,2] = chm[arrin[0,*,i], arrin[1,*,i]]
	endfor

	;right
	arrin = ulonarr(2,n,n_elements(seeds[0,*]))
	arrin[1,*,*] = rebin(seeds[1,*],n,numseeds)
	arrin[0,*,*] = rebin(seeds[0,*],n,numseeds) + rebin(lindgen(n),n,n_elements(seeds[0,*]))

	for i=0L, numseeds - 1 do begin
		chmA = chm
		chmA[where(chm lt chm[seeds[0,*], seeds[1,*]]) - hthresh]=max(chm)
		val = min(chm[arrin[0,*,i],arrin[1,*,i]], minsc)
		length[i] += minsc
		if (arg_present(rays)) then r[*,i,3] = chm[arrin[0,*,i], arrin[1,*,i]]
	endfor

	if (arg_present(rays)) then rays = r
	return, length[*] / 4.0
end